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Introduction 


This  project  dealt  with  computational  methods  for  inverse  problems  related 
to  light  propagation  through  the  earth’s  atmosphere.  In  this  section  we 
describe  three  important  applications: 

1.  Image  deblurring. 

2.  Estimation  of  phase  and  object  from  phase  diversity  data. 

3.  Estimation  of  the  volume  refractive  index  profile  of  the  atmosphere 
from  wavefront  sensor  data. 

Of  fundamental  importance  in  each  application  is  the  aperture-plane 
phase,  or  wavefront  aberration,  cf>(x,y).  This  can  be  interpreted  as  the  devi¬ 
ation  from  planarity  of  the  wavefront  of  light,  at  planar  incidence  at  the  top 
of  the  earth’s  atmosphere,  after  it  has  propagated  through  the  atmosphere 
[10].  From  the  phase,  one  can  compute  the  point  spread  function  (PSF)  for 
a  ground-based  telescope.  Assuming  no  imperfections  in  the  optics  of  the 
telescope  and  an  incoherent  light  source,  the  on-axis  PSF  takes  the  form 

s=\T{Ae^)\\  (1) 

where  T  denotes  the  two-dimensional  Fourier  transform,  i  =  \/— T,  and  A 
denotes  the  aperture  function.  For  an  ideal  telescope,  A  =  1  inside  the 
telescope  aperture,  and  A  =  0  outside  the  aperture. 

From  the  PSF,  one  can  formulate  a  simple  model  for  the  blurred,  noisy, 
pixelated  image  recorded  by  a  charge  coupled  device  (CCD)  array  attached 
to  a  telescope, 


dij  =  f  [  s(xi  -  x!,  yj  -  y')  f  (x,  y)  dx' dy  +%•.  (2) 

- - ' 

(S  f)  (xi,yj) 

Here  /  denotes  the  light  source,  or  object,  and  represents  stochastic  noise 
in  the  image  formation  process.  The  indices  i,  j  denote  position  in  the  pixel 
array. 

We  can  now  describe  several  important  inverse  problems  arising  in  atmo¬ 
spheric  optics.  Perhaps  the  simplest  is  image  deblurring,  i.e.,  the  estimation 
of  the  object  /  in  (2)  given  pixelated  image  data  d^.  Due  to  the  ill-posedness 
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of  inverse  of  the  underlying  convolution  integral  operator  S  defined  in  equa¬ 
tion  (2),  the  image  deblurring  problem  is  highly  ill-conditioned.  For  this 
reason,  regularization  (see  [11])  must  be  applied.  This  results  in  a  large- 
scale,  mildly  ill-conditioned  optimization  problem.  If  prior  information  like 
nonnegativity  of  the  object  is  incorporated,  the  problem  becomes  constrained 
and  much  more  difficult  to  solve  numerically. 

In  practice,  the  PSF  s  in  the  model  (2)  as  well  as  the  object  /  may  be 
unknown.  In  this  case,  one  may  use  phase  diversity  [7,  9]  to  simultaneously 
estimate  both  the  phase  <j>  and  the  object  /.  The  key  idea  is  to  collect  suitable 
additional  data.  Denote  the  data  in  model  (2)  by 

d  =  s[<f>]  *  /  +  77.  (3)  • 

Here  s  =  s[(f>]  denotes  the  dependence  of  the  PSF  on  the  phase  </>  in  equation 
(1),  and  *  denotes  convolution  product.  In  two-channel,  single  frame  phase 
diversity,  one  collects  additional  data 

d '  =  s[<f>  +  0\*  f  +  r{.  (4) 

Here  6  denotes  a  known  phase  aberration,  for  example,  aperture  plane  de¬ 
focus,  imposed  before  the  second  image  is  formed.  From  the  two  images  d 
and  d',  one  can  estimate  both  the  phase  4>  and  the  object  / .  From  (1)  the 
dependence  of  PSF  s  on  phase  (j>  is  nonlinear.  Hence,  the  phase  diversity 
inverse  problem  is  nonlinear.  It  is  also  ill-posed. 

When  imaging  over  relatively  wide  fields  of  view,  the  situation  is  even 
more  complicated  in  that  the  PSF  s  may  depend  on  the  viewing  direction  9 
(the  model  (l)-(2)  is  accurate  only  for  narrow  fields  of  view).  The  aperture 
plane  phase  in  the  direction  6  can  be  described  as  a  line  integral 

y,9)  —  f  n(x(z,9),y(z\9),z)dz.  (5) 

Jz= 0 

Here  n  =  n(x,y,z )  denotes  the  refractive  index  of  the  atmosphere,  and 
(x (z;9),y(z;9),z)  denotes  a  light  ray  path  from  the  top  of  the  atmosphere 
(z  =  H)  to  the  bottom  (z  =  0).  A  device  called  a  wavefront  sensor  measures 
the  gradient  of  <f>.  The  inverse  problem  is  to  estimate  the  volume  refractive 
index  profile  n(x,y,z )  from  wavefront  sensor  measurements  in  several  direc¬ 
tions  6.  This  is  a  limited  angle  tomography  problem,  and  is  ill-posed.  Prior 
information  can  make  this  problem  much  more  tractable.  The  atmospheric 
refractive  index  profile  is  highly  correlated  with  atmospheric  turbulence.  At¬ 
mospheric  turbulence  is  concentrated  in  a  few  thin  layers,  and  it  can  be 
accurately  modeled  using  Kolmogorov  statistics.  See  [8]. 
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2  Objectives 

The  broad  objective  of  this  project  was  the  development  of  efficient  com¬ 
putational  algorithms  to  solve  inverse  problems  that  arise  in  atmospheric 
optics.  We  focused  on  the  three  specific  application  problems  enumerated  in 
the  introduction. 

These  three  problems  share  several  features.  Each  is  ill-posed,  so  regular¬ 
ization  is  required  to  obtain  accurate  solutions.  Each  problem  has  solutions 
which  are  functions  of  two  or  three  spatial  variables.  When  discretized,  these 
yield  ill-conditioned  systems  with  very  large  numbers  of  unknowns.  In  the 
case  of  phase  diversity,  these  systems  are  nonlinear,  and  in  the  case  of  image 
deblurring,  the  incorporation  of  prior  knowledge  makes  these  systems  non- 
negatively  constrained.  In  summary,  we  needed  to  numerically  solve  large, 
ill-conditioned,  possibly  constrained,  and  possibly  nonlinear  systems  of  equa¬ 
tions. 

To  compute  solutions  to  these  systems,  there  are  several  fairly  general 
“tricks  of  the  trade”,  e.g.,  quasi-Newton  methods  for  nonlinear  systems,  pro¬ 
jected  Newton  methods  for  nonnegatively  constrained  systems,  and  the  con¬ 
jugate  gradient  (CG)  method  for  large  linear  systems.  To  gain  robustness 
and  efficiency,  we  were  forced  to  develop  very  special  preconditioners.  Pre¬ 
conditioners  can  be  viewed  as  transformations  that  take  advantage  of  special 
problem  structure.  Our  preconditioners  are  described  in  the  next  section. 


3  Major  Accomplishments  and  New  Findings 

Major  accomplishments  of  this  project  included 

•  The  development  of  efficient  nonnegatively  constrained  optimization 
algorithms  for  image  deblurring.  This  includes  a  new  preconditioner 
based  on  a  sparse  approximation  to  the  blurring  operator.  See  [1,  2] 
for  details. 

•  The  development  of  efficient  preconditioners  for  the  joint  phase  and 
object  estimation  problem  in  phase  diversity.  These  preconditioners 
were  based  on  the  Hessian  of  the  (quadratic)  regularization  terms.  See 
[4].  This  paper  also  contains  a  careful  numerical  study  and  comparison 
of  trust  region  vs.  limited  memory  BFGS  methods  for  the  numerical 
solution  to  optimization  problems  arising  in  phase  diversity  estimation. 
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Data  for  this  study  was  obtained  from  the  US  Air  Force  Maui  Space 
Surveillance  Complex  in  collaboration  with  Dr.  David  Tyler. 

•  The  development  of  multigrid  preconditioned  conjugate  gradient  schemes 
for  volume  refractive  index  (turbulence)  estimation.  See  [3,  5,  6].  These 
schemes  make  efficient  use  of  the  layered  structure  of  the  atmospheric 
turbulence  profiles.  This  layered  structure  gave  rise  to  block-structured 
matrices.  We  employed  a  block  analogue  of  symmetric  Gauss-Seidel  it¬ 
eration  as  our  multigrid  smoother. 

•  The  publication  of  a  research  monograph  entitled  “Computational  Meth¬ 
ods  for  Inverse  Problems”  [11].  This  publication  presents  both  the  gen¬ 
eral  theory  and  specific  algorithms  for  the  solution  of  inverse  problems. 


References 

[1]  J.  M.  Bardsley,  Constrained  Optimization  Techniques  for  Image  Recon¬ 
struction ,  Ph.D  Thesis,  Department  of  Mathematical  Sciences,  Montana 
State  University,  May  2002. 

[2]  J.  M.  Bardsley  and  C.R.  Vogel,  A  Nonnnegatively  Constrained  Con¬ 
vex  Programming  Method  for  Image  Reconstruction,  submitted  to  SIAM 
Journal  on  Scientific  Computing,  July  2002. 

[3]  B.  L.  Ellerbroek,  Luc  Gilles,  and  C.R.  Vogel,  A  Computationally  Effi¬ 
cient  Wavefront  Reconstructor  for  Simulation  of  Multi-Conjugate  Adap¬ 
tive  Optics  on  Giant  Telescopes,  Proc.  SPIE  4839-116,  Adaptive  Optics 
System  Technologies  II  (2002). 

[4]  Luc  Gilles,  C.R.  Vogel,  and  J.M.  Bardsley,  Computational  Methods  for 
a  Large-Scale  Inverse  Problem  Arising  in  Atmospheric  Optics,  Inverse 
Problems,  18  (2002),  pp.  237-252. 

[5]  Luc  Gilles,  C.R.  Vogel,  and  Brent  Ellerbroek,  A  Multigrid  Precondi¬ 
tioned  Conjugate  Gradient  Method  for  Large  Scale  Wavefront  Recon¬ 
struction,  Journal  of  the  Optical  Society  of  America  A,  19  (2002),  pp. 
1817-1822. 


6 


[6]  Luc  Gilles,  B.  L.  Ellerbroek,  and  C.R.  Vogel,  Layer- Oriented  Multigrid 
Wavefront  Reconstruction  Algorithms  for  Multi-Conjugate  Adaptive  Op¬ 
tics,  Proc.  SPIE  4839-118,  Adaptive  Optics  System  Technologies  II 
(2002). 

[7]  R.  A.  Gonsalves,  “Phase  retrieval  and  diversity  in  adaptive  optics” ,  Op¬ 
tical  Engineering,  21  (1982),  pp.  829-832. 

[8]  J.  W.  Hardy,  Adaptive  Optics  for  Astronomical  Telescopes,  Oxford 
Press,  1998. 

[9]  R.  Paxman,  T.  Schultz,  and  J.  Fineup,  “Joint  estimation  of  object  and 
aberrations  by  using  phase  diversity” ,  J.  Optical  Soc.  Amer.  A,  9  (1992), 
pp.  1072-1085. 

[10]  M.  C.  Roggemann  and  B.  Welsh,  Imaging  Through  Turbulence,  CRC 
Press,  Boca  Raton,  1996. 

[11]  C.  R.  Vogel,  Computational  Methods  for  Inverse  Problems,  SIAM,  2002. 

4  Personnel  Supported  by  Grant 

•  The  PI:  Curtis  R.  Vogel,  Professor  of  Mathematics,  Montana  State 
University. 

•  2  Graduate  Research  Assistants:  Johnathan  Bardsley  and  Scott  Hyde. 
Bardsley  completed  his  Ph.D.  in  Mathematics  in  May  2002  under  the 
direction  of  the  PI.  He  recently  accepted  a  postdoctoral  research  posi¬ 
tion  with  the  Statistical  and  Applied  Mathematical  Sciences  Institute 
(SAMSI)  in  North  Carolina.  Scott  Hyde  is  a  Ph.D.  student  in  Statis¬ 
tics,  also  at  Montana  State  University. 

5  Publications 

1.  C.  R.  Vogel,  Negative  results  for  multilevel  preconditioners  in  image 
deblurring,  in  Scale-Space  Theories  in  Computer  Vision,  Edited  by  M. 
Nielsen,  P.  Johansen,  O.F.  Olsen,  and  J.  Weickert,  Springer- Verlag, 
1999. 


7 


2.  C.  R.  Vogel,  A  limited  memory  BFGS  method  for  an  inverse  problem  in 
atmospheric  imaging,  in  Methods  and  Applications  of  Inversion,  Edited 
by  P.  C.  Hansen,  B.  H.  Jacobsen,  and  K.  Mosegaard,  Lecture  Notes  in 
Earth  Sciences,  92  (2000),  pp.  292-304,  Springer- Verlag. 

3.  M.  Hanke,  J.  Nagy,  and  C.  R.  Vogel,  A  Quasi-Newton  approach  to 
nonnegative  image  restoration,  Linear  Algebra  and  Its  Applications, 
316  (2000),  pp.  223-236. 

4.  Luc  Gilles  and  P.  Tran,  Optical  Switching  in  Nonlinear  Chiral  Dis¬ 
tributed  Bragg  Reflectors  with  Defect  Layers,  Journal  of  the  American 
Optical  Society  B,  19  (2002),  pp.  630-639. 

5.  Luc  Gilles,  C.R.  Vogel,  and  J.M.  Bardsley,  Computational  Methods  for 
a  Large-Scale  Inverse  Problem  Arising  in  Atmospheric  Optics,  Inverse 
Problems,  18  (2002),  pp.  237-252. 

6.  Luc  Gilles,  C.R.  Vogel,  and  Brent  Ellerbroek,  A  Multigrid  Precondi¬ 
tioned  Conjugate  Gradient  Method  for  Large  Scale  Wavefront  Recon¬ 
struction,  Journal  of  the  Optical  Society  of  America  A,  19  (2002),  pp. 
1817-1822. 

7.  J.  M.  Bardsley,  Constrained  Optimization  Techniques  for  Image  Recon¬ 
struction,  Ph.D  Thesis,  Department  of  Mathematical  Sciences,  Mon¬ 
tana  State  University,  May  2002. 

8.  C.  R.  Vogel,  Computational  Methods  for  Inverse  Problems,  SIAM, 
Philadelpha,  2002. 

9.  J.  M.  Bardsley  and  C.R.  Vogel,  A  Nonnnegatively  Constrained  Convex 
Programming  Method  for  Image  Reconstruction,  submitted  to  SIAM 
Journal  on  Scientific  Computing,  July  2002. 

10.  B.  L.  Ellerbroek,  Luc  Gilles,  and  C.R.  Vogel,  A  Computationally  Effi¬ 
cient  Wavefront  Reconstructor  for  Simulation  or  Multi- Conjugate  Adap¬ 
tive  Optics  on  Giant  Telescopes,  Proc.  SPIE  4839-116,  Adaptive  Op¬ 
tics  System  Technologies  II  (2002). 

11.  Luc  Gilles,  B.  L.  Ellerbroek,  and  C.R.  Vogel,  Layer- Oriented  Multi¬ 
grid  Wavefront  Reconstruction  Algorithms  for  Multi- Conjugate  Adap¬ 
tive  Optics,  Proc.  SPIE  4839-118,  Adaptive  Optics  System  Technolo¬ 
gies  II  (2002). 


8 


6  Presentations 

•  Sept.  26-30,  1999.  Santa  Clara,  California.  Annual  Meeting  of  the  Op¬ 
tical  Society  of  America.  Attended  and  presented  contributed  paper 
entitled  “Fast  algorithms  for  nonnegatively  constrained  image  deblur¬ 
ring”  . 

•  March  7-8,  2000.  Maui,  Hawaii.  Workshop  on  Post-Detection  Process¬ 
ing  and  Inverse  Problems  in  Ground-Based  Imaging  at  the  Maui  High 
Performance  Computing  Center.  Attended  and  presented  talk  entitled 
“Computational  techniques  for  inverse  problems” . 

•  June  28-July  1,  2000.  Lake  St.  Wolfgang,  Austria.  Conference  on 
Inverse  Problems.  Attended  and  presented  paper  “Incorporating  sta¬ 
tistical  information  in  the  regularized  solution  of  ill-posed  inverse  prob¬ 
lems”. 

•  July  10-12,  2000.  San  Juan,  Puerto  Rico.  SIAM  2000  Annual  Meeting. 
Attended  and  gave  2  minisymposium  talks:  (i)  “Inverse  problems  in 
atmospheric  imaging I:  Mathematical  models”;  and  (ii)  “Regularization 
parameter  selection  from  noise  statistics” . 

•  Oct.  16-20,  2000.  Minneapolis,  Minnesota.  Participated  in  IMA  Work¬ 
shop  ’’Image  Processing  and  Low  Level  Vision”. 

•  Dec.  5-7,  2000.  Maui,  Hawaii.  Workshop  on  Post-Detection  Process¬ 
ing  and  Inverse  Problems  in  Ground-Based  Imaging  at  the  Maui  High 
Performance  Computing  Center.  Attended  and  presented  talk  entitled 
“Computational  Methods  for  Nonnegatively  Constrained  Image  Recon¬ 
struction”.  Postdoc  Luc  Gilles  also  attended  and  presented  a  talk  on 
joint  work,  “Numerical  Issues  Arising  in  Phase  Diversity  Estimation”. 

•  Dec  14-16,  2000,  Hong  Kong.  Workshop  on  Mathematics  in  Image 
Processing.  Attended  and  presented  invited  talk  entitled  “Statistically 
Based  Regularization  Parameter  Selection  Schemes  for  Image  Deblur¬ 
ring”  . 

•  June  18-22,  2001.  Conference  on  Applied  Inverse  Problems  in  Monteca- 
tini,  Italy.  Attended  and  gave  invited  talk  entitled  “Statistically-Based 
Regularization  Parameter  Selection  Methods  for  Regularized  Ill-Posed 
Problems”. 


9 


•  Nov.  5-7,  2001.  Albuquerque,  New  Mexico.  Attended  the  Optical 
Society  of  America’s  Conference  on  Signal  Recovery  and  Synthesis  and 
presented  talk  “A  comparison  of  computational  methods  for  phase- 
diverse  object  and  phase  reconstruction” . 

•  Jan.  15-17,  2002.  Maui,  Hawaii.  Attended  AFOSR-sponsored  Work¬ 
shop  on  Advanced  Topics  in  Space  Situational  Awareness  and  presented 
talk  entitled  “Some  thoughts  on  nonnegatively  constrained  image  re¬ 
construction”  .  Postdoc  Luc  Gilles  also  attended  and  presented  a  talk 
on  joint  work,  “Fast  Algorithms  for  Multiconjugate  Adaptive  Optics” . 

•  July  10-12,  2002.  Philadelphia,  Pennsylvania.  Attended  SIAM  An¬ 
nual  Meeting,  organized  minisymposium  “Computational  Methods  for 
Large-Scale  Inverse  Problems”,  and  presented  talk  “Multiconjugate 
Adaptive  Optics”. 

•  Aug.  22-23,  2002.  Maui,  Hawaii.  Presented  short  course  “Multiconju¬ 
gate  Adaptive  Optics”  at  the  Maui  Scientific  Research  Center. 

•  Aug.  24-26,  2002.  Waikaloa,  Hawaii.  Attended  SPIE  Conference 
on  Astronomical  Telescopes  and  Instrumentation  and  presented  talk 
“Multigrid  wavefront  reconstruction  algorithms  for  extreme  and  multi¬ 
conjugate  adaptive  optics” . 

7  Collaborative  Research  and  Transactions  at 
US  Air  Force  Laboratories 

During  the  course  of  this  project,  the  PI  visited  Air  Force  research  facilities 
on  the  island  of  Maui,  Hawaii,  four  times.  Three  of  these  visits  were  to  at¬ 
tend  annual  workshops  and  meet  with  fellow  researchers  at  the  Maui  High 
Performance  Computing  Center  (MHPCC).  The  fourth  visit  was  to  present 
a  short  course  on  multiconjugate  adaptive  optics  at  the  Maui  Scientific  Re¬ 
search  Center  (MSRC),  which  is  next  door  to  the  MHPCC.  The  Pi’s  MSRC 
contact  was  Dr.  Stuart  Jeffries. 

The  PI  also  collaborated  with  Dr.  David  Tyler  in  the  analysis  of  phase 
diversity  data.  This  data  was  collected  by  Dr.  Tyler  at  the  US  Air  Force 
Maui  Space  Surveillance  Complex  on  Mt.  Halleakala.  Results  are  described 
in  publication  [5]  above. 
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8  Inventions  or  Patent  Disclosures 

None. 

9  Summary 

This  project  dealt  with  computational  methods  for  inverse  problems  related 
to  light  propagation  through  the  earth’s  atmosphere.  It  dealt  specifically 
with  the  following  applications:  (1)  image  deblurring;  (2)  the  simultaneous 
recovery  of  object  and  phase  from  phase  diversity  data;  and  (3)  The  estima¬ 
tion  of  the  volume  refractive  index  profile  of  the  atmosphere  from  wavefront 
sensor  data. 

Each  of  these  three  application  problems  was  ill-posed,  and  regularization 
was  required  for  their  accurate  solution.  In  each  case,  the  implementation  of 
regularization  gave  rise  to  optimization  problems  which  were  large  and  ill- 
conditioned.  In  addition,  the  image  deblurring  problem  was  nonnegatively 
constrained,  and  the  phase  diversity  problem  was  nonlinear.  These  properties 
made  the  optimization  problems  very  difficult  to  solve. 

Several  fairly  general  numerical  techniques  were  available  for  these  prob¬ 
lems.  These  included  trust  region  and  limited  memory  BFGS  nonlinear  opti¬ 
mization  methods  and  conjugate  gradient  iteration  to  solve  the  third  problem 
(volume  refractive  index  estimation)  and  to  solve  linear  subproblems  arising 
in  the  other  two  applications.  The  key  to  obtaining  robust,  efficient  solutions 
was  to  employ  an  approach  called  preconditioning.  The  PI  and  his  collab¬ 
orators  developed  very  effective  special-purpose  preconditioners  for  each  of 
the  three  specific  applications. 

The  PI  collaborated  with  several  AFOSR-sponsored  researchers  on  the 
island  of  Maui,  Hawaii  during  this  project.  The  PI  applied  the  computational 
techniques  that  he  developed  to  image  data  collected  at  the  US  Air  Force 
Maui  Space  Surveillance  Complex  by  Dr.  David  Tyler.  The  PI  also  presented 
a  short  course  on  the  topic  of  multiconjugate  adaptive  optics  at  the  Maui 
Scientific  Research  Center. 

A  total  of  9  peer-reviewed  scientific  journal  articles  were  prepared  under 
this  project.  A  PhD  thesis  was  also  written  by  a  student  directed  by  the  PI 
and  supported  under  this  project.  In  addition,  a  research  monograph  entitled 
“Computational  Methods  for  Inverse  Problems”  was  written  by  the  PI.  This 
monograph  was  recently  published  by  the  Society  for  Industrial  and  Applied 
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Mathematics. 

Preprints  and  reprints  of  papers  prepared  under  this  project  can  be  down¬ 
loaded  directly  from  the  Pi’s  web  site  at 

http : //www .math .montana . edu/~vogel/ 
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